Controlling crystal symmetries in phase-field crystal 
o ■ models 

(N 

^1 Kuo-An Wu^, Mathis Plapp^ and Peter W Voorhees^ '^ 

I ^ Department of Materials Science and Engineering, Northwestern University, 

Evanston, Illinois 60208, USA 
CS| , ^ Physique de la Matiere Condensee, Ecole Polytechnique, CNRS, 91128 Palaiseau, 

France 

^ Department of Engineering Sciences and Applied Mathematics, Northwestern 

C/3 I University, Evanston, Illinois 60208, USA 

E-mail: mathis .plappSpolytechnique . f r , kuoan-wuQnorthwestern. edu, 
p-voorhees@northwestern . edu 

■ Abstract. 

We investigate the possibility to control the symmetry of ordered states in phase- 
^ field crystal models by tuning nonlinear resonances. In two dimensions, we find that 

j-^ . a state of square symmetry as well as coexistence between squares and hexagons can 

O I be easily obtained. In contrast, it is delicate to obtain coexistence of squares and 

liquid. We develop a general method for constructing free energy functionals that 
exhibit solid-liquid coexistence with desired crystal symmetries. As an example, we 
develop a free energy functional for square-liquid coexistence in two dimensions. A 
systematic analysis for determining the parameters of the necessary nonlinear terms 
is provided. The implications of our findings for simulations of materials with simple 
' cubic symmetry are discussed. 

od 

§ ; PACS numbers: 68.08.-p,45.70.Qj, 47.54.-r, 81.16.Rf 

. 1. Introduction 

X 

' The phase-field crystal (PFC) model has rapidly gained popularity in recent years 
as a potentially powerful tool to simulate the evolution of materials with atomistic 
resolution on time scales that are orders of magnitude larger than those that can be 
attained by molecular dynamics simulations. While it was originally [H |2] inspired by 
phenomenological equations developed in the general context of pattern formation [3l H] , 
it has now been established |5l El [7] that it can be obtained as a simplification of classical 
density functional theory (DFT) of freezing [8l[9]. 

To make the PFC model useful for practical applications in materials science, it has 
to be established how the model parameters need to be chosen in order to reproduce the 
physical properties of a given material as closely as possible. This is currently a very 
active area of research, and several contributions of the present volume are dedicated 
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to this subject. One of the most fundamental properties of a material is of course its 
crystal structure. As detailed below, the original PFC model lU [2] contains only two 
parameters: a scaled global density ip and a scaled temperature e. The latter indicates 
the distance to the critical point of the PFC model, which is located at ip = e = 0. 
For small e, the PFC model exhibits, besides the unstructured liquid phase, periodic 
solutions of nematic (stripe) and hexagonal symmetry in two dimensions, and of bcc 
symmetry in three dimensions. Recently, it was found that for larger values of e, fee 
and hep structure also exist in three dimensions fiU\ [TT| [12] . However, the parameter e 
also controls several other properties of the model, such as the width of the solid-liquid 
interfaces and their interfacial free energy. Since, to match a given material, these 
quantities need to be adjusted independently, more degrees of freedom are needed in 
the model. 

It turns out that it is not straightforward to obtain periodic ground states with 
symmetries that differ from the "natural" ones. This is due to the physics which 
governs the formation of periodic states in the PFC model: the liquid state is unstable 
with respect to the formation of periodic solutions if their wavelength falls within a 
narrow band centred around a characteristic length scale. The periodic density patterns 
that correspond to crystalline phases are stabilised by nonlinear terms, which lead to 
the interaction (resonance) of density waves with different unstable wave vectors. The 
hexagon/bcc and stripe patterns arise from the simplest nonlinearities involving triadic 
and quartic resonances of wave vectors with equal modulus, respectively, which explains 
their ubiquity in nature [1] . This also naturally explains why new structures appear for 
larger values of e [TOl [HI [12]: since the range of unstable wavelengths increases with e, 
new resonances become possible. 

A way to control the selection of crystal structures is thus to modify the nonlinear 
resonances. Two ways to do this have been explored in the literature, both inspired by 
earlier work in the physics of pattern formation. The first idea is to make two bands of 
wavelengths unstable, and to use the ratio between the two characteristic wavelengths 
as an additional geometric parameter which allows to stabilise patterns of a certain 
symmetry [13]. Using this approach, PFC models with a ground state of fee symmetry 
were recently developed [HI [15]. While these models are robust and versatile, they 
have two potential drawbacks, namely (i) in the analogy with DFT, the two distinct 
length scales should arise from two distinct peaks in the liquid structure factor; for 
the simplest version of the model, a structure factor that is very different from those 
typically measured has to be used, and (ii) the equation of motion for the density field 
contains spatial derivatives of up to 10th order, which makes simulations in geometries 
that cannot easily be handled in reciprocal space very cumbersome. 

The second approach, which will be our main focus, is to add new nonlinear terms 
to the model which modify the strength of the resonances and can thus favour patterns of 
different symmetries. Here, we will restrict our investigations to two dimensions, where 
the obvious missing state is a square pattern. The question of how a square pattern can 
arise from a rotationally invariant homogeneous state has been extensively studied in 
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the context of convection patterns ^16J, because squares can be observed under certain 
experimental conditions. It has been demonstrated that simple additional nonlinearities 
can stabilise squares [171 HE], and that squares can coexist with hexagons [191 120] . 
These phenomenological equations have already been used to study the transition from 
hexagonal to square symmetry upon a change of the control parameters, in a spirit very 
close to the one of the PFC model, before that name was actually coined [211 122] ■ The 
geometric properties of grain boundaries between two domains of square symmetries 
have also been studied [23]. Finally, isolated patches of squares in coexistence with the 
unstructured state ( "oscillons" ) have also been found in similar equations [2^ [25] . 

Here, we investigate whether it is possible to use this approach to design a simple 
and robust PFC model with square symmetry. We find that it is indeed straightforward 
to obtain squares as well as coexistence of squares and hexagons following the recipes 
found in the literature. In contrast, square-liquid coexistence can only be obtained with 
a combination of several nonlinear terms and carefully tuned parameters. We develop 
a systematic approach to determine suitable parameter ranges that yield stable square- 
liquid interfaces, and confirm our analytic calculations by numerical simulations. We 
also evaluate the elastic coefficients of the square solid for the parameters used in the 
simulations. 

The remainder of the paper is organised as follows. In Section |2l we define the 
model and recall the calculation of nonlinear resonances and the resulting selection of 
symmetry. In Section |3] we calculate the free energies of the various periodic patterns 
in a generalised PFC model in two dimensions, with particular emphasis on square- 
liquid coexistence. Section HI presents the comparison between our analytic results and 
numerical calculations; the elastic coefficients are evaluated in Section |5l and Section [6] 
gives a brief conclusion. 

2. The model 

2.1. General considerations 

For the purpose of exposition, we start with the standard PFC model [H [2] , which has 
the dimensionless free energy functional 



where is a dimensionless particle density measured from a constant reference value, and 
e is a constant. This form of the free energy functional was originally proposed by Swift 
and Hohenberg as a phenomenological description of patterns that emerge in various 
hydrodynamical systems [3|. In this context, e is the distance from the bifurcation 
threshold at which the unstructured (quiescent) state becomes unstable. We have used 
dimensionless units in which the modulus of the characteristic wave number is unity, 
and energy and time have been scaled by appropriate quantities as detailed in [6]. The 
manner by which this free energy functional and its parameters can be related to classical 
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DFT is detailed in [51 El E] and need not be repeated here. 

In the original Swift-Hohenberg equation, ip is a linear combination of the deviations 
of the temperature and velocity fields from the quiescent state. It is treated as a non- 
conserved order parameter, and thus the system is free to evolve towards the global free 
energy minimum. In contrast, in the PFC model the field ip is & local particle density 
and hence a conserved quantity. This provides an additional control parameter, namely, 
the average density ip = J ipdf/ J dr. 

The free energy of periodic states can be evaluated using the one-mode 
approximation, in which the density field is given by a sum of density waves, 

^(f) =ijj + 5^{f) = + J2^K- exp(ii?j- -r) = ^jJ + exp{iKj ■ f), (2) 

i j 

where Kj are the principal reciprocal lattice vectors of the considered structure (with 
\Kj\ = q), and A is the amplitude of the density waves, where we have used the fact 
that for a homogeneous solid all density waves have the same amplitude, \Ag] = A. 
This ansatz is inserted into the free energy functional, and the result is integrated over 
one unit cell. In this procedure, all terms that contain oscillatory exponential factors 
integrate out, and only the products of exponentials in which two, three, or four of 
the Kj^s sum up to zero contribute to the final result. Thus, non-trivial contributions 
arise from the triadic and quartic resonances. Whereas the latter can generate stripes 
(or squares, see below), hexagonal and bcc ordered phases are stabilised by triadic 
resonances because the principal reciprocal lattice vectors are able to form "triads" 
(i.e., closed triangles, for example, (110), (lOl) and (Oil) for bcc lattices). The cubic 
terms that are needed for triadic resonances to occur are generated by the expansion of 
the ip^ term: {i/j + 61^)^ = (Stp)"^ + 4ip{6ip)^ + . . . . Therefore, they are absent for ip = 
and become increasingly important with increasing This explains the sequence 
of phases found in the phase diagram of the PFC model [2] with increasing rolls 
(stripes) — )■ hexagons — )■ liquid in two dimensions, and rolls — )• hexagons — bcc — )■ liquid 
in three dimensions. In addition, since the free energy expressed in terms of A contains 
terms in A"^, A^, and A'^, the bifurcation from liquid to hexagons can be transcritical, 
and a first-order transition from liquid to hexagons/bcc with a finite coexistence zone in 
ip is possible. In contrast, for crystals with square, simple-cubic and face-centred-cubic 
lattices, no triads can be formed, and the free energy contains only terms in and 
A^. This makes the transition from liquid to squares supercritical, and it was shown by 
weakly nonhnear analysis |6l HI] that solid-liquid coexistence is impossible under these 
conditions. 

From the above considerations, we can conclude that the model has to fulfil three 
requirements in order to obtain square-liquid coexistence. First, the quartic resonances 
need to favour squares rather than stripes, second the transition from squares to liquid 
must be subcritical to give rise to solid-liquid coexistence, and third the square state 
must have a lower free energy than hexagons even in the presence of triadic resonances 
that occur for ip 0. We will address these questions in the following, and we will 
develop a systematic approach to determine coefficients that favour periodic states with 



Controlling crystal symmetries in phase-field crystal models 
cubic symmetry. 



5 



2.2. Anisotropic terms 

Let us first discuss nonlinear terms that can be used to favour squares in two dimensions. 
To this end, it is useful to consider rhombi, that is, a density field composed just of two 
sets of density waves with wave vectors ±Ki and ±7^2, which have equal magnitude q 
and form an arbitrary angle 6 between them, 

64'{r) = A [exp(ii?i ■ f) + exp(ii?2 ■ + c.c.]. (3) 

Obviously, there are no triadic resonances. The quartic resonances generated by the ip'^ 
term in equation ([1]) give a result that is independent of 6. Thus, this term does not 
favour any particular symmetry. 

One possibility of an "anisotropic" term, used in the literature [iHl [191 [20] , is a term 
proportional to t/'^A^?/'^, where A = is the Laplace operator. Here, we consider a 
more general set of nonlinear terms of the form g2n ip'^A'^ip'^ in the free energy functional, 
where g2n are constant coefficients. Following the procedure outlined above, we obtain 

92nJ dfiP^A'^tP^ = g2nj df (tP^A'^iStpy + 4tlj^6^A''6tp + 2tlj6^A''{6^f 

+2V'(5V^)'A"(5V^ + (S^fA'^iSijY) . (4) 
Substituting equation ([3]) into equation yields 

^ J df^'A^r = {-irg2n + |i?2p") 

+2A^ (|2i?i|2" + |2i?2p'' + 4|i?i + i^sT" + 4|i?i - iTsT")] , (5) 

where V = J df is the volume of one unit cell. Together with the vector addition 
properties illustrated in Figure [H 

\Ki + K2\ = 2gcos- 

|i?i-i?2| =2gsin^, (6) 
equation ([5]) can be rewritten as 

^ J df V^2A>2 = (-l)"(^2ng'" {iQ^'A' + 22^+2^^) + 8A' 

/isotropic ~l~ /anisotropic- (7) 

where 

/isotropic = (-l)"^72ng'" {iQ^P' A' + 2^^+^ A') (8) 

is independent of the angle ^, but depends on ip., and the anisotropic part of the free 
energy that depends on the angle and thus on the symmetry of the pattern is 

/anisotropic = 8(-l)"^72ng'M^ fcOS^" ^ + siu^" ^^j . (9) 



'cos^- ^ + sin^" ^- 
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Figure 1. Terms of the form ■0^A"'0^ give rise to a 6'-dependent free energy through 
the above vector addition operations. 

For n = 1, the dependence on 6 in the anisotropic part vanishes because of the 
trigonometric identity cos^(6'/2) + sin^(^^/2) = 1. However, for n = 2, the anisotropic 
part of the free energy is a simple function of 6 which is minimised by = 7r/2 if 
> and hj 9 = if g4 < 0. Similar results are obtained for n = 3; however, the 
square symmetry is favoured when qq < 0. It is interesting to note that for n > 4, 
the anisotropic part of the free energy contains higher order harmonics, which suggests 
a possible means to favour arbitrary angles between density waves by an appropriate 
combination of such terms. 

Another type of nonlinearity used in the literature [T71 |2T1 [22] to favour squares 
is IV^/^I^. Indeed, repeating the above calculation for rhombi with a term of the form 
S4I V'?/'|*^/4 with S4 a constant, we find 

y Jdf^\Vtlj\^ = S4q^A\5 + cos^e) (10) 

which can again be split in an isotropic and an anisotropic part. The anisotropic part 
is minimised for = 7r/2 if S4 > 0. Note that, in this case, the isotropic part does not 
depend on ip. 

The above calculations are valid both in two and three dimensions. In three 
dimensions, an angle of 7r/2 between density waves corresponds to a simple cubic 
structure, which can therefore be obtained in a straightforward manner using this 
method. 

2.3. Subcritical bifurcation 

For crystal structures without triadic interactions, since to lowest order the free energy 
can be expressed in even powers of A, a subcritical bifurcation is required for solid-liquid 
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coexistence. The generic form of a free energy that exhibits a subcritical bifurcation is 

f = aA^ - PA'^ + jA^ (11) 

where /3 and 7 are positive. The parameter a controls the growth rate of small 
perturbations of the liquid state (linear stability); the liquid is stable for a > 0. The 
ip'^ term in the standard PFC free energy functional shown in equation ([1]) results in a 
^/'-dependent quadratic coefficient. Within the one- mode approximation, for this model 
a ~ — e + (1 — q'^y + Stp"^ for density waves with wave number q. 

All nonlinear terms discussed so far yield contributions up to order A^. To obtain 
the necessary contribution in A^, some higher order nonlinearity needs to be added. 
An obvious possibility would be just which was indeed used in 125] to model 
oscillons. The problem with this choice is that for ip it generates additional cubic 
and quintic terms in the free energy that tend to lower the free energy for hexagonal and 
bcc phases. As a consequence, this nonhnearity tends to generate the sequence squares 
— )■ hexagons — )■ liquid in the phase diagram; this was the case for all the cases that we 
have investigated. 

An alternative choice, which will be made in the following, is to add sejVV'l^/S to 
the free energy, which generates terms of order A^, but no additional terms depending 
on i/j. As will be shown below, it is indeed possible to generate a subcritical bifurcation 
that favours squares by using both S4I V'?/'|''/4 and SqIViP^/G with S4 < and Sq > 0. 
However, the resulting square-liquid equilibrium is metastable, whereas the state of 
lowest energy are hexagons. In order to be able to adjust the respective energies of 
these states, we must also include the nonlinear terms of the form tp'^A'^ip'^ and tp'^A^ip'^ 
discussed previously. The procedure how to choose appropriate coefficients for all these 
terms will be detailed in the next section. 



3. Selection of crystal symmetry 
3.1. Free energies 

The free energy functional considered in the following is 



J df ( I [-6 + (V^ + 1)2] ^ + ^ + ^^'A^r + f V^'A^^^^ 



+ (12) 

In two dimensions, the competing patterns are rolls, squares and hexagons. The 
corresponding density fields in the one-mode approximation of equation ([2]) are 

V'roii(^) =i! + 2A cos(gx), 
V'square(^) = ^ + 2A (cos qx + COS qy) , 

2 COS COS — hcosgxj, (13) 

where we recall that q is the magnitude of the principal reciprocal lattice vectors. The 
free energy is computed by substituting equation ( IT^ into equation (IT^ and integrating 
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Figure 2. Free energy as a function of for (54, sg) — (^2,1), (54,56) — (0,0) and 
e = 0.1. A common tangent line between the liquid and square free energy curves is 
drawn (thick dashed line) to illustrate the existence of a wide coexistence region. 



over one unit cell, which yields 

/roll = /l + (-e + (1 - q'f + (3 + 2g,q' - 2g,q'')^^) 

- i (-3 - 3s4g' - 16(74?' + 64^76?') + y S6g'A^ (14) 

/square = /i + 2 (-6 + (1 - q^f + (3 + 2(74?" " 2(76?')^^') A^ 

- (-9 - 5s4g' - 32(74?' + 96(76?') A' + ^s^q^A\ (15) 

/hex = /l + 3 (-e + (1 - q^f + (3 + 2(74?' - 2(76?')V^') A^ 
+ 12i; [l + g^q^-g,q^) A' 

/ 45 27 \ 

- (-y - y S4?' - 84(74?' + 264(76?'j + 91^6?' A^ (16) 

where / is the free energy density, / = F/V, and //, = (— e + l)ip'^/2 + Tp^ /A is the free 
energy density of the liquid. The values of A and ? are determined by minimising the 
above free energy with respect to A and ?. One important difference with respect to the 
standard PFC model should be noted. In the latter, q appears only in the gradient term. 
Therefore, ? and A decouple, and the value of ? that minimises the free energy is always 
equal to unity. In contrast, the new terms introduce nonlinear couplings between A and 
?, which implies that both quantities depend in a nontrivial way on the parameters. 
The values of ? and A which minimise / are found numerically for each structure using 
the Newton- Raphson method and are then substituted into equations (IT^ . (fT5|) and 
(fT6|) to obtain the free energy density for different patterns. 
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Let us first check whether square-hquid coexistence can be obtained using only 
the terms proportional to S4 and Sq. In figure |2] we plot the free energies of squares 
and liquid for (54, Se) = (—2,1), ((74, (/e) = (0,0), and e = 0.1. The common-tangent 
construction shows that there indeed exists a wide coexistence region as a result of the 
bifurcation being subcritical. However, the square-liquid coexistence is metastable for 
these parameters, since both the roll and the hexagonal phases have a lower energy, as 
shown in figure [2l Note that, since the cubic term which generates hexagons breaks the 
symmetry between A and —A, there are two distinct hexagon solutions that correspond 
to different values of A, which are labelled as "Hex-1" and "Hex-2" respectively in the 
graphs. To modify this ordering of the free energies, we need to include the terms 
proportional to and gQ. In the following, we give a systematic procedure to choose 
appropriate values for the coefficients. 



3.2. Rolls vs. squares 

We first consider the stability of rolls and squares. As discussed previously, choosing 
the coefficients of the anisotropic terms of the right sign favours squares. However, for 
fixed coefficients, the relative energies of the phases may still depend on the value of 
ip. Since the solid-liquid coexistence occurs near the density at which the solid and 
liquid free energy curves intersect, stable square-liquid coexistence is only possible if 
the intersection of the square and liquid free energy curves occurs at a larger value of 
lipl than the one for rolls and liquid. The following analysis focuses on determining the 
conditions for the nonlinear coefficients for which this is the case. 

The principal reciprocal lattice vectors of rolls and squares cannot form triads; 
thus, the free energy density of rolls and squares has the general form for a subcritical 
bifurcation given by equation ( ITTj) . 

^ = /L+«A2-/3A^+7A^ (17) 

where the subscript 5* (solid) denotes rolls or squares. The solid and liquid free 
energy curves intersect when the free energy difference of solid and liquid becomes zero, 
Sf = fs — fL = CiA"^ — f3A^ + = 0. In addition, the free energy of the solid must be 
minimised with respect to A, which yields {d6f /dA) = 2aA — 4/3 A^ + 6'yA^ = 0. These 
two relations yield 

Comparing equations (IT^ . ( ITB]) and (ITB]) to equation (ITTI) . we obtain 

a = N[-e + {l-q^f + T^^), (19) 

where 

r = 3 + 2g,q^ - 2g^q\ (20) 

and iV = 1,2 and 3 for rolls, squares and hexagons, respectively. Since the exponential 
growth rate of perturbations with wave number q is proportional to —a, in order to 
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obtain an unconditionally stable liquid at large values of {ipl, a must have a finite lower 
bound which leads to 

r > 0. (21) 

The intersection of solid and liquid free energy curves is obtained by combining 
equations ( fTSj) and ( !T9|) . which yields 

The wave numbers q that minimise the free energy of rolls and squares are different. 
However, for large values of Sg this difference is small so that we can safely neglect it. 
Under this assumption, the intersection point of solid and liquid free energy curves is 
solely dictated by the value of /3^/(4A^7). The condition that |'?/'squarel > iV'Sul then 
yields inequalities between the coefficients which are given in the appendix. The most 
important result is that this condition can only be satisfied if Qq is non-zero. 



3.3. Hexagons 

For hexagons, the free energy contains a cubic term due to triadic interactions, 

/hex = fL + aA^ + tA^ - f3A' + 7A^ (23) 



where 



a = 3(-e+(l-g2)2 + r^2 



T = 6ip{T - 1) 
45 2 T 

{3= - — - —s,q^ - Mg^q" + 2Ug,q^ 

7 = Slseg' (24) 

In the standard PFC model, the cubic term changes the bifurcation from supercritical to 
transcritical which not only make hexagons the favoured phase, but also makes hexagon- 
liquid coexistence possible in the limit of a weakly first-order freezing transition. The 
cubic term plays a similar role in the subcritical case: it lowers the free energy for 
hexagons by an amount proportional to In order to make square symmetries 

favourable, we require the cubic coefficient r to be small for if) close to the square- 
liquid coexistence region. Then, the subcritical bifurcation analysis developed in the 
previous section still gives a good estimate of the intersection points for the free energy 
curves of hexagons and liquid and can thus be used to evaluate the relative stability of 
squares and hexagons. 

This requirement on r and the condition for a stable liquid are illustrated 
graphically in figure |3]for ((74, (^e) = (~3.0, —2.1). The condition that the liquid remains 
linearly stable for large values of requires that F > 0. The cubic coefficient r is 
a function of q and proportional to the distance between the solid and dashed lines 
in figure [3l We set (54, Sq) = (—25, 600) and e = 0.001 so that these parameters 
satisfy the conditions for a subcritical bifurcation and for iV^gquarel > l^roiil shown in 
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Figure 3. T = 3 + 2g4,q^ - 2^69^ (solid line) for (54,56) = (-3,-2.1). The condition 
for a stable liquid is satisfied since F is positive. The cubic coefficient r = 6V'(r — 1) 
is proportional to the distance between the solid and dashed lines. 



equations ( IA.5I) and ( IA.7I) in the appendix. The wave number q is about 0.98 near the 
coexistence region, for which the cubic coefficient r is small as shown in figure |3l We can 
then estimate iV'hexl using the above analysis for a subcritical bifurcation by assuming 
r 0; details can be found in the appendix. For the parameters chosen above, both 
and l^/'hexl have smaller values than l^/'squarel- These parameters are used in the 
following numerical simulations. 



4. Comparison of analytical solutions and numerical simulations 

For a conserved order parameter ip, the evolution equation for the free energy functional 
shown in equation ( 1T2|) is 



dt Sij 

= ([-e + (V^ + 1)2]^ + + g^^A""^^ + ^?6M'^' 

- S4V ■ iWl^V^) - seV ■ (I WI'VV')) . (25) 

Since the main focus of this paper is the stability of crystal symmetries at equilibrium, 
we use a nonlocal globally conserved dynamics [26j to accelerate the search of equilibrium 
solutions. 



6F 1 



df 



SF 



(26) 



dtp 

'dt^ ' 6iIj ' V J dtp' 
We set the simulation parameters in two dimensions to be {g4,,ge) = (—3.0, —2.1), 
(s4, Se) = (—25,600), e = 0.001, and the grid spacings Ax and Ay close to 0.5. 
Numerical simulations are carried out using equation (12^ for rolls, squares and hexagons 



Controlling cry; 



12 



0,04 



0,03' 



|A| 0,02 
0,01 



1 ' 

• * 


1 






^ A A 
♦ 


^ f o 
♦ ♦ 


-- Rolls 

— Squares 
Hex-1 
Hex-2 

1 


/ 

/' 

/ 

■ 





-0,08 -0,06 ^0,04 -0,02 

Figure 4. Comparison of ttie magnitude of A that minimises the free energy as a 
function of t/" for different patterns. The analytical solutions are plotted as lines. The 
numerical simulations for rolls, squares, Hex-1 and Hex-2 are plotted as squares, circles, 
diamonds and triangles, respectively. 



in one unit cell with initial conditions listed in equation f lT^ . The free energy is evaluated 
numerically using equation (fT2|) after the numerical solution reaches a steady state. Since 
the wave number q that minimises the free energy depends on t/), the simulation for each 
pattern at a fixed value of is repeated with a different system size until the minimum 
of the free energy is found. 

A comparison of the amplitudes of the principal reciprocal lattice vectors obtained 
from the numerical simulations and the corresponding analytical solution is shown 
in figure HI In the numerical simulations, the amplitudes of the principal reciprocal 
lattice vectors are computed using the Fourier transform. For rolls, the one-mode 
approximation and the numerical simulations are in good agreement. This is reasonable 
because the next-nearest reciprocal lattice vectors of rolls are (20), which is far from the 

(10) principal reciprocal lattice vectors. This makes them difficult to excite, and the one- 
mode approximation is quite accurate. In contrast, for squares and hexagons, the higher- 
order reciprocal lattice vectors are closer to the principal reciprocal lattice vectors (e.g., 

(11) for squares), and thus the one-mode approximation and the numerical simulations 
do not agree as well. Thus, an analytical calculation including higher-order modes is 
required to give more accurate predictions. Nevertheless, the one-mode approximation 
gives a good qualitative prediction of the relative values of the free energies at the end 
point of the solid solution branches, as well as for the intersection points of solid and 
liquid free energy curves. 

The free energy difference f — Jl is plotted in figure [51 This difference is so small 
that it would be difficult to distinguish the free energy curves if they were plotted as 
in figure [2j The one-mode approximation and numerical simulations show the same 
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Figure 5. Comparison of tlie free energy densities as a function of -0 for different 
patterns. Symbols and lines as in figure ID 



ordering of the intersection points of solid and liquid free energy curves, as in figure HJ 
namely l^/'squarel > l^hexl > I Varolii- addition, for the parameters chosen above, rolls 
and hexagons are metastable states and the square lattice is the ground state for a wide 
range of \ip\ until it loses its stability to liquid at \ip\ = 0.073. To illustrate the influence 
of triadic interactions on the free energy calculation, a three-mode approximation that 
considers (10), (11) and (12) reciprocal lattice vectors for the square lattice is computed 
and shown in figure [5l 

To simulate solid-liquid coexistence, we set the simulation parameters to Ax = 
Ay = 0.53 on a system of size = 384Ax and Ly = 192A?/ with periodic boundary 
conditions. The initial condition for the density is 

^ = l{^+ tanh (x - + " ^'"''^ ~ tO ^'^'^^ 

where ips is the one-mode approximation of square lattices as shown in equation ( !T3|) . 
and ipL is the constant density of the liquid. The average densities of solid and 
liquid are chosen to be close to ^square' which is determined numerically. The initial 
amplitude of the square lattice is set to the value obtained from the steady state 
square lattice simulations at the same average density. Uniformly distributed random 
fluctuations of the magnitude of half the square pattern amplitude are applied initially 
to examine the stability of the square pattern. The simulations show that the square- 
liquid coexistence is stable against the initial random fluctuations, and the equilibrium 
solid-liquid coexistence is shown in figure El Furthermore, it can be seen in figure |6]that 
the square-liquid interface displays different spatial decay rates for the different density 
waves as predicted by the classical density functional theory of freezing. In particular. 
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Figure 6. Square- liquid coexistence simulated using the free energy functional shown 
in equation ([T^. 



the density wave of (10) decays more slowly into the liquid than the density wave of (01) 
for the {10} interface. The directional dependence of the spatial decay rate of density 
waves into the liquid was shown to be a main determinant of the anisotropy of interface 
properties, as discussed in [6| [27]. 



5. Elastic constants 



We briefly investigate here the elastic properties of our model. This is important because 
the analysis of the two-mode model presented in [14j predicts that the shear modulus 
for square patterns is zero in the limit of small e if only the mode corresponding to the 
principal reciprocal lattice vectors (10) is active. Only the addition of the second mode 
corresponding to the second reciprocal lattice vectors (11) ensures a finite and positive 
shear modulus. Since our model is based on a single unstable mode, it is important to 
evaluate the shear modulus. 

To determine the elastic constants, we follow the lines of earlier work [2], IHj and 
consider perturbations of the density field given by equation (fT3|l . For shear, bulk and 
deviatoric deformations, we use 

^shear = + 2A (cos q{x + ^y) + COS qy) , (28) 



Anik = tp + 2A f cos -j^^r^ + cos -^Y^p^ j ' (29) 



and 



V'deviatoric = ^ + 2A [cOiS. J^^^^ + {1-^) ) ' ^^^^ 

respectively, where ^ is the strain. These expressions are then inserted in the free 
energy functional ( !T2|) . which can be explicitly evaluated in terms of g. A, and ^. We 
set the wave number q equal to the wave number of the reference (unstrained) state. 
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and determine A by minimising the free energy of the strained state for each value of 
^ numerically using again the Newton-Raphson method. This value is then used to 
evaluate the free energy density, as a function of ^. The difference with the free energy 
of the unstrained state, A/, can be related to the elastic constants. To lowest order in 
C,, we have 

A f - 

^ J shear 2 

A/bulk = {Cu + c,2)e 

A/deviatoric = (Cll — ^12)^^- (31) 

We have evaluated the free energies of the strained states for the same parameters 
as the numerical simulations shown previously, except that we set ip = for simplicity. 
The results indeed display the expected quadratic behaviour for ^ ^ 1, and the elastic 
constants are, in our dimensionless units, 

Cu = 7.8 X 10"^ 
Cu = 1.1 X 10-^ 

C44 = 1.63 X 10^^ (32) 

The shear modulus is /i = C44 = 1.63 x 10~^, and the (two-dimensional) bulk modulus 
B = (Cii + Ci2)/2 = 3.85 X 10~^. The ratio of the bulk modulus and the shear modulus 
is about 24. 

Thus, for the parameters chosen here the square patterns generated by our model 
are rather "soft" with respect to shear deformation, but have a perfectly well-defined 
finite shear modulus. Several remarks can help to understand this fact. A shear 
deformation corresponds, in reciprocal space, to a change of the angle between reciprocal 
lattice vectors, whereas their length remains unchanged to first order in ^. Since we 
have found that the free energy depends on this angle if anisotropic terms are present 
according to equation (|9]), this directly gives a finite contribution to the shear modulus. 
Note, however, that this is not the only contribution. Indeed, since the change in the 
angle modifies the strength of the nonlinearities which saturate the density waves, the 
amplitude A also depends on the deformation, which modifies the contributions of all 
the other terms in the functional. Furthermore, the nature of the bifurcation also comes 
into play. For a subcritical bifurcation, A can be appreciable even for small e, which 
means that the anisotropic terms, though proportional to or A^, may be comparable 
to the quadratic contribution which determines the elastic constants in the standard 
PFC model. 

Note that the above calculation of the elastic constants uses the one-mode 
approximation, which is not very accurate in view of the results shown in figure HI 
However, it is sufficient for our main purpose here, which was to demonstrate that 
the shear modulus is finite. For a more quantitative evaluation of the elastic constants, 
multi-mode calculations or direct numerical computations are mandatory. Furthermore, 
we have evaluated the elastic constants here only for one specific set of coefficients. How 
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they depend on the various parameters in the free energy functional is an interesting 
subject for further investigations. 

6. Conclusions 

We have presented a general method of constructing a free energy functional for the 
PFC model that exhibits solid-liquid coexistence with desired crystal symmetries. We 
have demonstrated that the crystal symmetries in the PFC model can be controlled 
through additional nonlinear terms. In particular, we have examined the influence of 
nonlinearities such as 5'2n'^^'^"'^^ and S2m|VV^p'" on crystal symmetries. Besides the 
squares that have been investigated here, other structures can potentially be obtained 
by using terms with n > 4, since these terms contain higher order harmonics that can 
be used to favour well-defined angles between density waves. 

We have also presented a systematic procedure for constructing free energy 
functionals that exhibit square-liquid coexistence as the ground state. This requires 
(i) a subcritical bifurcation from the liquid to the square state, and (ii) a suitable 
combination of the nonlinear coefficients which favours squares over rolls and, at the 
same time, makes triadic resonances small enough to avoid hexagons. Our analytical 
results obtained in the one-mode approximation are borne out qualitatively by numerical 
simulations: in both simulations and analytics the square lattice is the ground state until 
it loses its stability to the hquid at large {ipl. The one-mode approximation, however, 
does not predict accurately the amplitude and free energy for square lattices. A multi- 
mode approximation is needed to quantitatively describe the square symmetry since 
the principal lattice vectors of the square lattice are strongly coupled to higher order 
modes. Nevertheless, the one-mode analysis provides a good guideline for determining 
coefficients of nonlinearities to obtain the desired crystal symmetries. 

We have obtained stable square-liquid coexistence, and the solid exhibits a finite 
shear modulus, which makes this model suitable for simulations. A drawback is that 
several nonlinear terms are needed, and that a delicate balance between the coefficients 
needs to be respected. We have only tested one specific set of coefficients, but we expect 
that the properties of the model (elastic constants, interfacial properties) can be varied 
by exploring the parameter space of the various coefficients. It should also be noted 
that our survey of potential nonlinear terms is by no means complete. However, the 
relation between the coefficients and the resulting model properties is highly non-trivial, 
which implies that the process of finding a set of coefficients that yield a match with 
some desired materials properties is likely to be cumbersome. Furthermore, our final 
equation of motion ( l25|) is of 8th order in space, which is only two orders lower than the 
one of the two-mode model |14] . 

A highly interesting perspective is the extension of this work to three dimensions. 
The calculation of the anisotropy performed in section |2] remains valid in three 
dimensions, which means that it is straightforward to obtain simple cubic structures. 
Furthermore, since it has been found that the free energies of bcc, fee, and hep phases 
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in three dimensions in the standard PFC model are not very different for a certain range 
of e |TOl [m [12] , the addition of anisotropic terms may be used to control their relative 
stability. 
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Appendix A. Inequalities for the coefficients in the free energy functional 

Here, we give the details concerning the calculation of the intersection points between 
liquid and solid free energies and the resulting conditions for the coefficients in the free 
energy functional. We start with rolls and squares. As outlined in section 3.2, we need to 
compare the values of (4A^7), where /3, and 7 are the coefficients of equation ( |T7I) for 
a generic subcritical bifurcation. Comparing equations (IT^ and (ITS!) to equation (fT7|) . 
we obtain the coefficients /3 and 7 for rolls and squares, 

/^rou = ^(-3 - 354?"^ + AiS), 

7roll = y■56?^ (A.l) 
/^square — ( 9 5^45' + ^) , 

56 

Tsquare = y'^S'?^, (A. 2) 

where 

E = - 32(74?' + 96(76?', 

X^ = i-lQg,q^ + QAg,q^)/E. (A.3) 

The subcritical bifurcation requires the quartic coefficient to be negative, which 
according to equation (ITTi) implies 

Aquarc > 0. (A.4) 

This determines an upper bound for S4, 

s^q' < (A.5) 
5 

The condition that the intersection of the square and liquid free energy curves occurs 
at a higher value of \ip\ than that of rolls and liquid (i.e., |'?/'squarel > iV^roiil); yields 

A^square ^ /^roll (A 6) 



S^square 47roll 

which can be reduced to the inequality 

99 - 25S + 42AiH - Jnl < s^q^ < 99 - 25H + 42AiH + Jvl^, (A.7) 
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where 

= 70(5AiH - 3S + 12)2. g^. 

The solution of s^q"^ exists only if the upper bound of S4q'^ obtained from equation ( lA.Sp 
is greater than the lower bound obtained from equation (]A.7p . which yields 

(-9 + S)/5 > 99 - 25S + 42AiS - VtO |5AiH - 3S + 12|. (A.9) 

The above inequality holds if (a) Ai < 3/5 and S > 12/(3 — 5Ai) or (b) Ai > 3/5 and 
S < 12/(3 — 5Ai). It is helpful to rewrite the expression of Ai as 

A, = 5 + 1 . !m1±^. (A.10) 

For the case that = 0, we have Ai = 3/5 — 1/10 < 3/5 and the solution of S4 only 
exists if S > 12/ (3 — 5Ai) > 0. However, this requires g^ to be negative which contradicts 
the condition for stable liquid shown in equation fl2T]) . Thus it is essential to include the 
nonlinear term ip'^A^ip'^ in the free energy functional so that the free energy functional 
exhibits the desired properties of (i) stable liquid, (ii) subcritical bifurcation, and (iii) 

IV^squarel > I Varolii • 

For hexagons, the coefficients a, /3, and 7 are given by equation The coefficient 
r of the cubic term is assumed to be small, such that this term can be neglected. Then, 
the condition that iV'squarel > iV^hexl requires 



/^square Phex 



> J^^^ (A.ll) 



SO'square 1 27he 

yielding the inequality 

- 45 + 65S + 24A2S - ^2 < SAq"" < -45 + 65S + 24A2S + ^fi^, (A.12) 

where 

A2 = (84^4?' - 264(76g')/H, (A.13) 

and 

^2 = — (IOA2S + 27S-I8)'. (A.14) 
9 

Together with equation f[0|) . we find that either (a) A2 > -27/10 and S > 18/(27 + 
I5A2) or (b) A2 < -27/10 and S < 18/(27 + I5A2) has to be fulfilled in order to make 

hex I ■ 
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